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The algebraic Bethe Ansatz is a prosperous and well-established method for solving one- 
dimensional quantum models exactly. The solution of the complex eigenvalue problem is thereby 
reduced to the solution of a set of algebraic equations. Whereas the spectrum is usually obtained di- 
rectly, the eigenstates are available only in terms of complex mathematical expressions. This makes 
it very hard in general to extract properties from the states, like, for example, correlation functions. 
In our work, we apply the tools of Tensor Network States to describe the eigenstates approximately 
as Matrix Product States. From the Matrix Product State expression, we then obtain observables 
like correlation functions directly. 



The Coordinate Bethe Ansatz as originally in- 
vented by Bethe to give an exact solution to the one- 
dimensional antiferromagnetic Heisenberg model, re- 
duces the complex problem of diagonalizing the Hamil- 
tonian to finding the solutions of a set of algebraic equa- 
tions. Once solutions to these algebraic equations are 
found - numerical approaches to find them efficiently ex- 
ist in many cases - the eigenvalues are known exactly. 
However, the eigenstates are available only as a complex 
mathematical expressions. This makes it insuperable to 
get interesting properties out of the states - like their 
entanglement characteristics or their correlations. 

A complementary approach is the Algebraic Bethe 
Ansatz [2[ (also known as "inverse scattering method"). 
In this approach, the scattering matrix (i?-tensor) is in 
focus. Based on this matrix, the Hamiltonian is derived 
and eigenstates are constructed. The "inverse" problem 
consists in finding the scattering matrix that represents 
favored Hamiltonian. As in case of the Coordinate Bethe 
Ansatz, the eigenvalue problem is reduced to solving a 
set of algebraic equations. From the solutions, the eigen- 
values are obtained directly. The eigenstates are avail- 
able only in terms of a non-contractable tensor network. 
This makes exact calculations of expectation values in 
general unfeasable. However, because of their inherent 
structure, it has been proven to be possible to calculate 
the norm and the scalar product between Bethe states 
exactly [HQ. 

In this paper, we take advantage of the fact that the 
Bethe-eigenstates have the form of a tensor-network 0, 
\^ . The calculation of correlation functions then requires 
the contraction of a tensor network such as the one de- 
picted in Fig. [TJ We contract the tensor network approx- 
imately using a similar method as for time evolution in 
the Density Matrix Renormalization Group (DMRG) 
0]. Finally, we end up with a Matrix Product State 
(MPS) [1, Hoi - Hi) from which we can extract expectation 
values like correlation functions directly. We show for 
the case of the antiferromagnetic Heisenberg model and 
the XXZ model with both periodic and open boundary 
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FIG. 1. Tensor network describing correlations {(Jx^cfx^) 
with respect to a Bethe-eigenstate. 



conditions that correlations can be obtained for 50 sites 
with good precision. 

The Bethe-eigenstates for the antiferromagnetic 
Heisenberg model and the XXZ model with periodic 
boundary conditions are obtained as products of oper- 
ators B(fjLj) applied on a certain vacuum state | vac), i.e. 

|*(^i,.,.,/xm)> = B(fi 1 ) • • • B(fi M )\vac). (1) 

The parameters {fij} are thereby solutions of Bethe equa- 
tions and the £?(/ij)'s play the role of creation operators 
of down-spins. The vacuum corresponds to the state with 
all spins up. B(/jl) has the form of the Matrix Product 
Operator (MPO) Q 
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with k,l e {0,1}, of = I k)(l\ (0 =t, 1 =i) and £f(/i) 
being 2x2 matrices dependent on the parameter \i. The 
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product of operators B(/a\) • • • B(/jlm) can be read as the 
contraction of the set of 4-index tensors [Cf {/Aj)] r r , with 
respect to a rectangular grid, as shown in Fig. El Thereby, 
r, r', k and I label the left, right, up and down- indices, 
respectively. The non-zero entries of the tensor [Cf {/Aj)\ r r i 
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[£J( M )]? = b((i) 



with b(/i) = 1/(1 + /x), c(/i) = /i/(l + /x) for the Heisen- 
berg model and 6(/i) = sinh(2^)/ sinh(/i + 2iry), c(/i) = 
sinh(/i)/(sinh(/i-h2i77) for the XXZ-model (77 is related to 
the inhomogenity A in the XXZ model via A = cos(2t?)). 

The calculation of expectation values with respect 
to such a Bethe-eigenstate is a considerably com- 
plex problem, because it requires the contraction 
of the tensor network depicted in Fig. [U This 
tensor network represents the correlation function 

(^(/ii, . . . ,/xm) WaPvx } | ^(/ii, • • • ,M )• It is com- 
posed of the network for | . . . , /am) ) shown in 
Fig. [2j the vertically mirrored network with conjugated 
tensor entries representing ( ^(/ii, . . . , /am) I and two op- 
erators cfx^ and cr^ squeezed in-between at sites i and j. 
A tensor-network with such a structure also appears in 
connection with the calculation of partition functions of 
two-dimensional classical systems and one-dimensional 
quantum systems [9] and the calculation of expectation 
values with resp ect to Projected Entangled Pair States 
(PEPS) [l5l4l7jl. The complexity to contract this net- 
work scales exponentially with the number of rows M or 
columns N (depending on the direction of contraction), 
which render practical calculations infeasible. 

To circumvent this problem, we attempt to perform 
the contraction in an approximative numerical way: the 
main idea is to consider the network in Fig. [2] as the 
time-evolution of the state with all spins up by the 
evolution operators B(/ai), . . . , B(/am) to the final state 
I ^(/ii, . . . , /am) )• After each evolution step, the state re- 
mains a MPS, but the virtual dimension is increased by 
a factor of 2. Thus, we approximate the MPS after each 
evolution step by a simpler MPS with smaller virtual 
dimension. Of course, caution has to be used, because 
the operators B(/aj) are not unitary and the intermedi- 
ate states of the evolution can be non-physical (i.e. they 
might have to be represented by MPS with high virtual 
dimension). The Bethe- network, however, bears an ad- 
ditional structure which can be taken advantage of: all 
evolution operators B(/aj) commute, such that the oper- 
ators B(/Aj) can be arbitrarily ordered. This allows us to 
choose the optimal ordering for the evolution with inter- 
mediate states that are least entangled. We will discuss 
this in detail the following. 

The algorithm consists of M steps m = 1, . . . , M: in 
the first step, the vacuum- state | vac ) is multiplied by the 
MPO to form the initial MPS | *i ). In step m > 



1, the product | \I/ m ) = B(/A m ) \ ^ m -i ) is approximated 
by the MPS | \I/ m ) that has maximal bond-dimension D 
and is closest to | \I/ m ). In other words, we try to solve 
the minimization problem 



K:= \* m )-\* m ) 



Min. 
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by optimizing over all matrices of the MPS |^ m ). 
This minimization problem also appears in the context 
of numerical calculation of expectation values with re- 
spect to PEPS [l5l-[l7|. the calculation of partition func- 
tions [sl and (imaginary) time-evolution of ID quantum 
systems [H . We discuss it in detail in the supplementary 
part. In this way, the MPS-approximation of the Bethe- 
state is obtained for m = M. Thereby, {/ii, . . . , /am} 
are the solutions of the Bethe-equations. The error of 
the approximation is well-controlled in the sense that 
the expectation- value of the energy can always be calcu- 
lated with respect to the approximated MPS | ^m ) and 
compared to the exact energy available from the Bethe- 
ansatz. 

In order to make the algorithm more efficient, we 
take into account the "creation operator" -property of the 
MPO B(/i): Since each MPO B{/a) creates one down- 
spin, the MPS I \£ m ) at step m contains exactly m down- 
spins. Explicitly, the MPS reads 

l*m)= (0\(0\A kl ---A kN \0)\m)\k 1 ,...,k N ) 

ki • • -fcjv 

with matrices A k being block-diagonal in the sense that 
(a\(s\A k \P)\s') = [A k }f s ,. a and /? are the virtual 
indices that range from to D — 1 (with D being the vir- 
tual dimension of the state). One the other hand, s and 
s' are the symmetry indices that transfer the information 
about the number of down-spins from left to right. The 
local constraint that guarantees this information trans- 
fer is s' = s + k. This constraint determines the blocks 
[w4 fe ]l^/ that are non-zero and allows a sparse storage 
of the state. The left boundary-state ( | and the right 
boundary-state | m ) fix the total number of down-spins 
of the MPS to m. The optimization problem [2] at step m 
then consists in approximating the state | \I/ m ) with m 
down-spins by a state | \I/ m ) that also has m down-spins. 
This leads to a gain of a factor of m in time and memory. 

Furthermore, there is a (mathematical) degree of free- 
dom that can be used to improve the approximation. 
This degree of freedom is due to the commutativity- 
property of B(/a): since [B(/a), B{y)\ — for all /a 
and z/, the ordering of the B(/ij) f s in the Ansatz (pp) 
is completely arbitrary. This is relevant insofar as the 
entanglement-properties of the intermediate states | \I/ m ) 
(1 < m < M) are concerned. The intermediate states 
are a priori no physical ground states, i.e. there is no 
reason for them to lie in the set of MPS with low bond- 
dimension. However, as we see numerically, there is al- 
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FIG. 2. Tensor network constituting the Bet he eigenstate of FIG. 4. Tensor network constituting the Bet he eigenstate 
the Heisenberg model or XXZ model with periodic boundary of the Heisenberg model or XXZ model with open boundary 



conditions. 



conditions. 
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FIG. 3. Half-chain entropy as a function of the evolu- 
tion step for the ground state of the N = 12 Heisenberg 
model with periodic boundary conditions. The lower line 
is obtained with the orderings (435261), (435216), (432561), 
(432516), (345261), (345216), (342561), (342516); the up- 
per line is obtained with the orderings (531624), (426153), 
(351624), (246153). 

ways an ordering such that the intermediate states con- 
tain as little entanglement as possible. This ordering we 
then use for calculating the approximation. That the 
ordering has a formidable effect can be gathered from 
Fig. [3] for the example of the ground state of the 12-site 
Heisenberg antiferromagnet with periodic boundary con- 
ditions. Here, the half-chain entropy of | \I/ m ) is plotted 
as a function of m for the best and the worst ordering. As 
it can be seen, the entropy is highest at the intermediate 
steps and decreases at m —> M when the state becomes 



a physical ground state. 

Up to now, we have considered the Bethe Ansatz for 
the XXZ model with periodic boundary conditions (and 
models in this class). In case of open boundary condi- 
tions, the Bethe Ansatz has the same form as in (pQ), 
merely the creation Operators are not single MPOs, but 
products of two MPOs |, [3, [l9|: 

l 

s=0 

Bi- 3 (fjL) has the property to create 1 — s down-spins, 
whereas B s (fi) creates s down-spins (s G {0,1}), such 
that B(/i) is a creation operator for exactly one down- 
spin, as before. The tensor- network representation for 
the Bethe-state with open boundary conditions is shown 
in Fig. 2J It contains twice as many rows as the tensor- 
network for periodic boundary conditions, which makes 
the contraction more challenging, in principle. However, 
as we see numerically, after a multiplication with a MPO- 
pair B(/i), the Schmidt-rank of the state increases only 
by a factor of 2 - not 4, as expected. Thus, the numerical 
effort to contract the tensor-network is similar. 

Using the previously described method, we have ob- 
tained results for the Heisenberg model with periodic 
boundary conditions and the XXZ-model with open 
boundary conditions. 

In case of the Heisenberg model, we have investigated 
the two-spinon excited states with total spin 5 = 1 and 
total z-spin S z = 1. The energies of these states obtained 
by the Bethe- Ansatz as a function of the momentum for 
N = 50 spins are plotted in the right inset in Fig. With 
respect to these states, we have calculated the correlation 




FIG. 5. Structure factor for the ground state (dotted line) 
and three selected two-spinon excited states of the periodic 
N = 50-Heisenberg chain as a function of the wave-vector q 
obtained from approximated Bethe-MPS with D — 1000. The 
three states are marked in the inset a the rhs which shows all 
energies and momenta of the two-spinon excited states with 
total spin S = 1 and total z-spin S z = 1 The inset at the lhs 
depicts the relative error in the energy as a function of D for 
the three selected two-spinon states. 




FIG. 6. Structure factor S(ir) for selected two-spinon excited 
states of the N = 50 Heisenberg antiferromagnet as a function 
of the momentum and the excitation energy. The used cutoff 
is D = 1000. 

functions {(J r z (J s z ) and the corresponding structure factor 

The structure factor S z (q) as a function of q for the three 
selected excited states (that are marked in the right in- 
set) can be gathered from Fig. In order to judge the 
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FIG. 7. Structure factor S(q) of the ground states of the 
N = 50 XXZ-model with open boundary conditions for se- 
lected values of A obtained from approximated Bethe-MPS 
with D — 1000. The dashed black line marks the free fermion 
case A = and the dotted black line indicates the Heisenberg 
limit A = 1. The inset shows the overlap with MPS of virtual 
dimension D = 400 obtained from DMRG calculations. 

accuracy of the approximated Bethe-MPS obtained with 
our method, we have compared the the expectation val- 
ues of the Hamiltonian with respect to these MPS to the 
energies obtained from the Bethe-Ansatz. For the case 
N = 50, the so obtained relative error plotted as a func- 
tion of D for the selected excited states can be gathered 
from the left inset in Fig. [5] As can be seen, the er- 
ror decreases with increasing D, but also with decreasing 
excitation energy. In Fig. [6j the structure factor at the 
point q = 7r, i.e. the squared staggered magnetization, 
is plotted as a function of the excitation energy and the 
momentum. Evidently, the excited states of the lowest 
branch show the highest staggered magnetization. 

In case of the XXZ-model with open boundary con- 
ditions, we have studied the correlations of the ground 
state. To get an impression of the quality of the ob- 
tained result, we computed the overlap with MPS states 
obtained from DMRG calculations. For the case of 50 
spins, the overlap as a function of D can be gathered 
from the inset of Fig. [7] for different values of A. The 
structure factor S(q) as a function of the wave- vector q 
for different values of A evaluated with respect to the 
approximated Bethe-MPS with D = 1000 is plotted in 
the main part of Fig. [71 The structure factor obtained 
for this D only deviates marginally from the structure 
factor obtained from the DMRG calculation. 

Summing up, we have presented a method for approx- 
imative calculation of correlation functions with Bethe- 
eigenstates. For this, we make use of the fact that a Bethe 
eigenstate is a product of MPOs applied to a MPS. We 
systematically reduce the virtual dimension after each 



multiplication and obtain an MPS with small virtual di- 
mension that can be used for the calculation of any ex- 
pectation value. We have shown the effective operation 
of this method by applying it to the Heisenberg anti- 
ferromagnet with periodic boundary conditions and the 
XXZ model with open boundary conditions. We have 
obtained results for the structure factor of ground- and 
excited states and compared our ground state-results to 
DMRG calculations. 
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Supplementary Material 

In this supplementary material, we describe the numerical method used in our article in detail. 

Numerical State Approximation 

The main building block of the algorithm is the approximation of MPS | \I/ m ) with a fixed number of m down-spins 
and virtual dimension D by a MPS | \I/ m ) that also has m down-spins and virtual dimension D < D in a way that 
the distance (j2j) between the two states is minimal. The state | \I/ m ) reads 

l*m>= E (0|(0|^---^|0)|m)|fc 1 ,...,fc N > 

ki---k]\f 

with (a\(s \A k | f3 ) | s' ) = [^J^J, and s' = s + k. The MPS | \I/ m ) has the same structure, but may be inhomogeneous, 
i.e. the matrices may be site-dependent: 

l*m)= £ (0K0|^ 1 ---i^|0)|m)|fc 1 ,... ) fc N ) 

ki---kjsr 

We obtain the starting point for the optimization by setting | \I/ m ) equal to | \I/ m ) and performing Schmidt- 
decompositions successively for bonds (1, 2) to (N — 1, TV), where we keep at each bond only the D largest Schmidt- 
coefficients. Due to the conservation of the number of down-spins, the Schmidt-decomposition with respect to one 
bond (j, j + 1) can be written as 

m r s 

s=0 7 =0 

Thereby, the states {| )\j = 1, . . . , T s } are states acting on the left block (sites 1, . . . , j) with a fixed number of s 
down-spins. In addition, these states satisfy the orthonormality constraint ( <1> 7 | <1> 7 , ) = 5 77 /. In an analogous manner, 
|| q>rn-s ^ _ i 5 . . . 5 Y s } are orthonormal states acting on the right block (sites j + 1, . . . , TV) with a fixed number 
of m — s down-spins. The singular values corresponding to the partitioning of s down-spins in the left block and m — s 
down-spins in the right block are {<j^|7 = 1, . . . , T s }. 

The Schmidt-decomposition can be obtained by transforming the MPS into the form 

l*m>= £ (0\(0\- ■■ AfXAfX ■■■\0}\m}\k 1 ,...,k N ), 

ki • • -kjsf 

with Ai fulfilling the local constraints Y^k(^i)^i = 1 for i = 1, . . . , j and ^2k ^H^iY = 1 for i = j + 1, . . . , iV. 
The local constraints guarantee the orthonormality of the states in the left and right block. E is a diagonal matrix 
containing the Schmidt-coefficients. The transformation into this form is always possible due to the gauge invariance 
of MPS (I . 

For i < j, the way to meet the local constraints is by QR-decompositions of the matrices A\ defined as [Af ]p sk = 
[jl- 1 ]^/, successively for i = 1, ... ,j: Af = Qf Rf . The matrix Qf is an appropriate replacement for Af , because 
the orthogonality of Qf makes the local constraint satisfied. To keep the MPS invariant, the block-diagonal matrix Ri 
defined as = [Rf }pdss' must be used to update A k +1 to RiA k +1 . Only at the last step i = j, the update must 

be omitted. For i > j, LQ-decompositions are performed in an analogous manner: successively for i = TV, . . . , j + 1, 
the matrices [»A|]^ s / fc = 1*4?]^/ are decomposed as LfQf. The matrices Qf are then used to replace Af and Lf update 
A k _ 1 to A k _ 1 L i (with [Li]p s s , = [Lf]p5 3S /). As before, the update must be omitted at the last step i = j + 1. 

The Schmidt-decomposition is now obtained by performing a singular value decomposition of the product of the 
"left-over" update-matrices R? and Lj +1 \ RjLj +1 = U s T i s V s . U s and V s are unitary, as well as their block-diagonal 
extensions U and V defined as [?7]^|, = [U s ]p5 ss > and [V]^, = [U s ]p5 ss >. Because of their unitarity, they can be used 
to update Aj to AjU and Aj +1 to VA k +1 without spoiling the orthonormality constraint. The matrix E s contains 
the Schmidt-coefficients {<r^|7 = 1, . . . ,T S } related to the partitioning of s down-spins in the left block and m — s 
down-spins in the right block. With the definition of E as = [S S ]^(5 SS /, we have obtained the desired form of 

the MPS. 
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A sensible way to reduce the dimension of bond (j, j + 1) is keeping only the largest D Schmidt-coefficients and 
setting all others to zero. Thus, by defining a projector P, such that PEP^ is the D x Z)-matrix containing the D 
largest Schmidt-coefficients and updating Aj to A k P^ and Aj +1 to PE^ +1 , we obtain the favoured MPS with 
reduced bond-dimension. 

Performing the Schmidt-decompositions with successive projections for all bonds (1,2) to (N — l,iV) gives a MPS 
that is a fairly good starting point for the optimization problem (j2j). A further improvement is possible by optimizing 
the quantity K locally, i.e. by optimizing K with respect to the matrices {Aj \ k = 0, 1} at one site j, and keeping all 
other matrices constant. This has already been discussed extensively in [9||8[. The main idea is that K is a quadratic 
function of Aj, such that it can be written as 



K = const + Ys^yNfA)' + 



kk' k 

The minimum with respect to Aj is achieved for those values of Aj solving the system of linear equations 

k> 

The matrix Mf k ' is a function of A*, i ^ j, and it is equal to the identity if the constraints ^2 k (A k y A k = t for 
i = - 1 and E/c A fc (A fe ) f = 1 for i = j -h l,...,iV. are fulfilled. These constraints, however, can always 

be imposed, as argumented earlier. By performing the local minimization for j sweeping between 1 and N until 
convergence of K, the global minimum of K is (usually) obtained and we have found the optimal approximation 
| ty m ) with maximal bond-dimension D to the state | \I/ m ). 



